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Abstract 

An efficient computational algorithm to price financial derivatives is presented. It is 
based on a path integral formulation of the pricing problem. It is shown how the path 
integral approach can be worked out in order to obtain fast and accurate predictions 
for the value of a large class of options, including those with path-dependent and 
early exercise features. As examples, the application of the method to European and 
American options in the Black-Scholes model is illustrated. A particularly simple 
and fast semi-analytical approximation for the price of American options is derived. 
The results of the algorithm are compared with those obtained with the standard 
procedures known in the literature and found to be in good agreement. 
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1 Introduction 



The classical theory of option pricing is based on the results found in 1973 by 
Black and Scholes [1] and, independently, Merton [2]. Their pioneering work 
starts from the basic assumption that the asset prices follow the dynamics 
of a particular stochastic process (geometric Brownian motion), so that they 
have a lognormal distribution [3,4]. In the case of an efficient market with no 
arbitrage possibilities, no dividends and constant volatilities, they found that 
the price of each financial derivative is ruled by an ordinary partial differen- 
tial equation, known as the Black-Scholes-Merton (BSM) formula. In the most 
simple case of a so-called European option, the BSM equation can be explic- 
itly solved to obtain an analytical formula for the price of the option [3,4]. 
When we consider other financial derivatives, which are commonly traded in 
real markets and allow anticipated exercise and/or depend on the history of 
the underlying asset, the BSM formula fails to give an analytical result. Ap- 
propriate numerical procedures have been developed in the literature to price 
exotic financial derivatives with path-dependent features, as discussed in de- 
tail in Refs. [3,5,6]. The aim of this work is to provide a contribution to the 
problem of efficient option pricing in financial analysis, showing how it is pos- 
sible to use path integral methods to develop a fast and precise algorithm for 
the evaluation of option prices. 

The path integral method, which traces back to the original work of Wiener 
and Kac in stochastic calculus [7,8] and of Feynman in quantum mechanics [9], 
is today widely employed in chemistry and physics, and very recently in finance 
too [10-14], because it gives the possibility of applying powerful analytical and 
numerical techniques [15]. Following recent studies on the application of the 
path integral approach to the financial market as appeared in the econophysics 
literature (see Refs. [12,14] for a comprehensive list of references), this paper 
is devoted to present an original, efficient path integral algorithm to price 
financial derivatives, including those with path-dependent and early exercise 
features, and to compare the results with those obtained with the standard 
procedures known in the literature. 

The paper is organized as follows. In Section 2 the basic ideas of the classical 
theory of option pricing are summarized, discussing the computational com- 
plexity associated to the evaluation of the price of a path-dependent option 
and reviewing the standard numerical procedures adopted in the literature. In 
Section 3 the path integral approach to option pricing is described and ana- 
lytically developed in order to obtain an efficient procedure for the calculation 
of the transition probability associated to a given stochastic model of asset 
evolution. Theoretical and computational details to obtain fast predictions 
for path-dependent options are also described. As applications of the method, 
numerical results for European and American options in the BSM model are 
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given in Section 4, together with comparisons with results known in the lit- 
erature. A particularly simple and very quick semi-analytical approximation 
for the price of an American option is derived in Section 5, by exploiting 
the possibility of anticipated exercise for any time before the expiration date. 
Conclusions and possible perspectives are drawn in Section 6. 



2 Option pricing: theory and numerical procedures 

2.1 Classical theory and path- dependent options 

The basic ingredient for the development of a theory of option pricing is a 
suitable model for the time evolution of the asset prices. The assumption 
of the BSM model is that the price S of an asset is driven by a geometric 
Brownian motion and verifies the stochastic differential equation (SDE) [3,4] 

dS = fiSdt + aSdw, (1) 

which, by means of the ltd lemma, can be cast in the form of an arithmetic 
Brownian motion for the logarithm of 5* 

d(\nS) = Adt + adw, (2) 

where a is the volatility, A = (// — cr 2 /2), /z is the drift parameter and w is 
the realization of a Wiener process. By virtue of the properties of a Wiener 
process, eq. (2) may be written as 

d(ln S) = Adt + aeVdt, (3) 

where e follows from a standardized normal distribution with mean and 
variance 1. 

Thus, in terms of the logarithms of the asset prices z' = In S',z = InS, the 
conditional transition probability p(z'\z) to have at the time t' a price S' under 
the hypothesis that the price was S at the time t < t' is given by [4,10][3] 

[z'-{z + A{t'-t))f \ 

2a 2 (t'-t) \' 1 1 



1 The correct way to indicate conditional transition probabilities is p{z' , t'\z, t). We 
omit the times in order to simplify the notation. 



p{z'\z) 



2n(t' - tW 



exp 
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which is a gaussian distribution with mean z + A(t' — t) and variance a 2 (t' — t). 
If we require the options to be exercised only at specific times t i: i = 1, ■ ■ ■ ,n, 
the asset price, between two consequent times t^_i and t iy will follow eq. (3) 
and the related transition probability will be 



with At = ti — ti-i. 

A time-evolution model for the asset price is strictly necessary in a theory of 
option pricing because the fair price at time t = of an option O, without 
possibility of anticipated exercise before the expiration date or maturity T (a 
so-called European option), is given by the scaled expectation value [3] 

O(0) = e- rT E[0(T)}, (6) 



where r is the risk-free interest and E[-] indicates the mean value, which can 
be computed only if a model for the asset underlying the option is understood. 
For example, the value O of an European call option at the maturity T will be 
maxjS'r — X, 0}, where X is the strike price, while for an European put option 
the value O at the maturity will be max{X — St, 0}. It is worth emphasizing, 
for what follows, that the case of an European option is particularly simple, 
since in such a situation the price of the option can be evaluated by means of 
analytical formulae, which are obtained by solving the BSM partial differential 
equation with the appropriate boundary conditions [3,4]. On the other hand, 
many further kinds of options are present in the financial markets, such as 
American options (options which can be exercised at any time up to the ex- 
piration date) and exotic options [3], i.e. derivatives with complicated payoffs 
or whose value depend on the whole time evolution of the underlying asset 
and not just on its value at the end. For such options with path-dependent 
and early exercise features no exact solutions are available and pricing them 
correctly is a great challenge. 

Actually, in the case of options with possibility of anticipated exercise before 
the expiration date, the above discussion needs to be generalized, by introduc- 
ing a slicing of the time interval T. Let us consider, for definiteness, the case 
of an option which can be exercised within the maturity but only at the times 
t\ = At, t% = 2 At, . . . , t n = nAt = T. At each time slice £j_i the value 
of the option will be the maximum between its expectation value at the time 
ti scaled with e~ rAt and its value in the case of anticipated exercise Oj_ x 0. 
If Si-i denotes the price of the underlying asset at the time we can thus 

2 For example, the value O of a call option in the case of exercise at the time 
will be max{Si_i — X, 0}, X being the strike price. 
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write for each % = 1, . . . , n 



0;_i(<Si-i) = max 



{or^Si-^e^EPilSi-r]}, 



(7) 



where E'fOjlS'j-i] is the conditional expectation value of 0j, i.e. its expectation 
value under the hypothesis of having the price Si-i at the time In this 
way, to get the actual price Co, it is necessary to proceed backward in time 
and calculate O n -i, . . . ,Oi, where the value O n of the option at maturity is 
nothing but 0%(S n ). It is therefore clear that evaluating the price of an option 
with early exercise features means to simulate the evolution of the underlying 
asset price (to obtain the Oj) and to calculate a (usually large) number of 
expectation conditional probabilities. 



2.2 Standard numerical procedures 

To value derivatives when analytical formulae are not available, appropriate 
numerical techniques have to be advocated. They involve the use of Monte 
Carlo (MC) simulation, binomial trees (and their improvements) and finite 
difference methods [3,5]. 

A natural way to simulate price paths is to discretize eq. (3) as 



which is correct for any At > 0, even if finite. Given the spot price So, i- e - the 
price of the asset at time t — 0, one can extract from a standardized normal 
distribution a value e^, k — 1, . . . , n for the random variable e to simulate one 
possible path followed by the price by means of eq. (8): 



Iterating the procedure m times, one can simulate m price paths {(So, S[ J , S-f', 
. . . , S^p = S^) : j = 1, . . . , m}, to which apply the procedure exemplified in 
Section 2.1 and evaluate the price of the option. In such a MC simulation of 



In S(t + At) - In S(t) = AAt + aeVAl, 



or, equivalently, 



S(t + At) = S(t) exp \AAt + ae^/At] , 



(8) 



S(kAt) = S{{k - l)Af) exp \AAt + ae k VAt . 
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the stochastic dynamics of asset price (Monte Carlo random walk) the mean 
values E[Oi\Si-i\, i — 1, . . . , n can be simply obtained as 

+ o {2) + ■■■ + o {m) 

m 



with no need to calculate transition probabilities because, through the extrac- 
tion of the possible e values, the paths are automatically weighted according 
to the probability distribution function of eq. (5). 

Unfortunately, this method leads to an estimated value whose numerical error 
is proportional to vrT 1 ! 2 . Thus, even if it is powerful because of the possibil- 
ity to control the paths and to impose additional constrains (as it is usually 
required by exotic and path-dependent options), the MC random walk is ex- 
tremely time consuming when precise predictions are required and appropriate 
variance reduction procedures have to be used to save CPU time [3]. 
This difficulty can be overcome by means of the method of the binomial trees 
and its extensions (see [3] and references therein), whose main idea stands in a 
deterministic choice of the possible paths to limit the number of intermediate 
points. At each time step the price Si is assumed to have only two choices: 
increase to the value uSi,u > 1 or decrease to dSi,0 < d < 1, where the 
parameters u and d are given in terms of a and At in such a way to give the 
correct values for the mean and variance of stock price changes over the time 
interval At. Also finite difference methods are known in the literature [3] as an 
alternative to time-consuming MC simulations. They provide the value of the 
derivative by solving the differential equation satisfied by the derivative, by 
converting it into a difference equation. Although tree approaches and finite 
difference methods are known to be faster than the MC random walk, they 
are difficult to apply when a detailed control of the history of the derivative 
is required and are also computationally time consuming when a number of 
stochastic variables is involved [3]. It follows that the development of efficient 
and fast computational algorithms to price financial derivatives is still a key 
issue in financial analysis. 



3 The path integral method 



The path integral method is an integral formulation of the dynamics of a 
stochastic process. It is a suitable framework for the calculation of the transi- 
tion probabilities associated to a given stochastic process, which is seen as the 
convolution of an infinite sequence of infinitesimal short-time steps [10,15]. 
For the problem of option pricing, the path integral method can be employed 
for the explicit calculation of the expectation values of the quantities of finan- 



6 



cial interest, given by integrals of the form [10] 

EpilSi-t] = J dzfflizilzi-JOiie*), (9) 



where z — In 5" and p(zi\zi-i) is the transition probability. E[Oi\Si-i] is the 
conditional expectation value of some functional Oi of the stochastic process. 
For example, for an European call option at the maturity T the quantity 
of interest will be max{SV — X, 0}, X being the strike price. As already 
emphasized, and discussed in the literature [3,5,6,11,14], the computational 
complexity associated to this calculation is generally great: in the case of 
exotic options, with path-dependent and early exercise features, integrals of 
the type (9) can not be analytically solved. As a consequence, we demand 
two things from a path integral framework: a very quick way to estimate the 
transition probability associated to a stochastic process (3) and a clever choice 
of the integration points with which evaluate the integrals (9). In particular, 
our aim is to develop an efficient calculation of the probability distribution 
without losing information on the path followed by the asset price during its 
time evolution. 



3.1 Transition probability 

The probability distribution function related to a SDE verifies the so-called 
Chapman-Kolmogorov equation [4], i.e. the relation 

p(z"\z') = J dzp{z"\z)p{z\z'), (10) 



which states that the probability (density) of a transition from the value z' (at 
time t') to the value z" (at time t") is the "summation" over all the possible 
intermediate values z of the probability of separate and consequent transitions 
z' — > z, z — > z". 

As a consequence, if we consider a finite time interval [f, t"] and we apply a 
time slicing, by considering n + 1 subintervals of length At = (t" — t')/n + 1, 
we can write, by iteration of eq. (10) 

dz n p(z"\z n )p(z n \z n - 1 ) ■ ■ -p(zi\z'), 



+oo +oo 



p(z"\z')= J--- Jc 



dz\ 



— CO — oo 
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which, thanks to eq.(4), can be written as 

+00 +00 



J- ■ ■ ldz x ■■■dz n 



1 



—00 —00 



yJ(2na*At) 



n+1 



• cxp 



In the limit n — > 00, At — > such that (n+l) At = (t"—t r ) (infinite sequence of 
infinitesimal time steps), the expression (11), as explicitly shown in Ref. [10], 
exhibits a Lagrangian structure and it is possible to express the transition 
probability in the path integral formalism as a convolution of the form [10] 



1 

p(z",t"\z',t') = J P[cr- 1 z]exp< - J L(z(r),k(T);r)dT 



where L is the Lagrangian 



L{~z{r)/i{r)-r) = 



2a 2 



z(t) - A 



and the integral is performed (with functional measure T>[-}) over the paths 
z(-) belonging to C, i.e. all the continuous functions with constrains z(t') = z', 
z(t") = z". As carefully discussed in Ref. [10], a path integral is well defined 
only if both a continuous formal expression and a discretization rule are given. 
As done in many applications, the Ito prescription is adopted in the present 
work. 

A first, naive evaluation of the transition probability (11) can be performed 
via Monte Carlo simulation, by writing eq. (11) as 



+00 +00 



P (A t> V , f) =/.../ ft Hh-j^ «p {-5^25 V ~ <* + ^)] 2 } .(12) 



—00 —00 



in terms of the variables g, L defined by the relation 

d9k = v^m exp [Zk ~ {Zk ~ i + AAt)]2 }' (13) 

and extracting each gi from a gaussian distribution of mean Zk-i + AAt and 
variance a 2 At. However, as we will see, this method requires a large number of 
calls to obtain a good precision. This is due to the fact that each gi is related 
to the previous ft-i, so that this implementation of the path integral approach 
can be seen to be equivalent to a naive MC simulation of random walks, with 
no variance reduction. 
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By means of appropriate manipulations [15] of the integrand entering eq. (11), 
it is possible, as shown in the following, to obtain a path integral expression 
which will contain a factorized integral with a constant kernel and a consequent 
variance reduction. We will refer to this second implementation of the method 
as path integral with importance sampling. 

If we define z" = z n+ \ and y k = z k — kAAt, k — 1, . . . , n, we can express the 
transition probability distribution as 

T T 1 f 1 n+l } 

■■■ dy 1 ---dy n — r • exp \ -— — J2 [Vk - Vk-i? \, (14) 



in order to get rid of the contribution of the drift parameter. Now let us extract 
from the argument of the exponential function a quadratic form 



n+l 



h/k ~ Vk-i? = 2/o- 2 2/i2/o + y\ + vl ~ 2yiy 2 + • • • + vl+i 



k=i 



y l My + [yl - 2y l y + y 2 n+l - 2y n y n+1 ], 



(15) 



by introducing the n-dimensional array y and the nxn matrix M defined as 

(l -10 \ 

-12-10 •••0 



y = 



2/2 



M = 



-1 2 -1 ••• 

•••-12-10 

-12-1 

\0 -12 



(16) 



where M is a real, symmetric, non singular and tridiagonal matrix. In terms 
of the eigenvalues mi of the matrix M, the contribution in eq. (15) can be 
written as 



y l My = w'O* MOw = w'M d w = ^m^ 2 , 



(17) 



by introducing the orthogonal matrix O which diagonalizes M, with w; 
OijUj. Because of the orthogonality of O, the Jacobian 



J = det 



dwi 
dy k 



det|Cy , 
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of the transformation y^ — > Wk equals 1, so that EHLi dwi = 117=1 dy%- Thanks 
to eqs. (16)-(17), and after some algebra, eq. (15) can be written as 



n+l 



k=l 



i=l 



n 

i=i 



W.; 



(VoOli + Vn+lOni) 



1 2 



mi 



+ y 2 + y 2 n+i -± iy ° Oll+ l n+lOm)2 .m 



i=l 



mi 



Now, if we introduce new variables hi obeying the relation 



dhi = \ I — exp 



2na 2 At 



2a 2 At 



W,; 



(yoOu + y n +iO r 



ITLi 



>dw i: (19) 



it is possible to express the finite-time probability distribution p(z"\z') as 

+oo +oo 



/ / fl d V*-n=^ 



exp , 

v /(2-^A/)"- : I 2n T - 

1 



— oo — oo 



-oo +oo 



1 n+1 ) 

k=l J 



11 dWi 



-oo — oo 



1=1 



y / (27ra 2 At)"+ 1 



e -(s/d+<+i)/2- 2 Ai. 



exp 



2a 2 At 



+ oo -|-oc 

/ /II 



dhi 



rrii Wi 



1 



(VoOli + Vn+lOni) \ (y O U + Vn+lOnif 



mi 



mi 



— oo — oo 



J i=i ,/27rcr 2 Atdet(M) 



• exp 



2a 2 At 



yn+lO n i) 



i=l 



mi 



-.(20) 



Equation (20) is one of the main results of the present work. Actually, the 
probability distribution function, as given by eq. (20), is an integral whose 
kernel is a constant function (with respect to the integration variables) and 
which can be factorized into the n integrals 



+oo 

/ 



dhi exp 



1 (yoOu + y n +iO 

r, 



2a 2 At 



rrii 



(21) 



given in terms of the hi, which are gaussian variables that can be extracted 
from a normal distribution with mean (yoOu + y n +iO n i) 2 /m^ and variance 
a 2 At/rrii. Differently to the first, naive implementation of the path integral, 
now each hi is no longer dependent on the previous hi-i, and importance 
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sampling over the paths is automatically accounted for. 

It is worth noticing that, by means of the extraction of the random variables 
hi, we are creating price paths, since at each intermediate time ij the asset 
price is given by 



Therefore, the path integral algorithm can be easily adapted to the cases in 
which the derivative to be valued has, in the time interval [0, T] , additional 
constraints, as in the case of interesting path-dependent options, such as Asian 
and barrier options [3]. 

The results of the two realizations of the path integral method here discussed 
will be compared in Section 4. 



3.2 Integration points 

Thanks to the method illustrated in Section 3.1, a powerful and fast tool to 
compute the transition probability in the path integral framework is available 
and it can be employed if we need to value a generic option with maturity 
T and with possibility of anticipated exercise at times U = iAt (nAt = T). 
As a consequence of this time slicing, one must numerically evaluate n — 1 
mean values of the type (9), in order to check at any time tj, and for any value 
of the stock price, whether early exercise is more convenient with respect to 
holding the option for a future time. To keep under control the computational 
complexity and the time of execution, it is mandatory to limit as far as possible 
the number of points for the integral evaluation. This means that we would 
like to have a linear growth of the number of integration points with the time. 
Let us suppose to evaluate each mean value 



with p integration points, i.e. considering only p fixed values for Zi. To this 
end, we can create a grid of possible prices, according to the dynamics of the 
stochastic process as given by eq. (3) 



n 



Si = exp {^2 O ik h k + iAAt}. 



(22) 



fc=i 




z(t + At) - z(t) = In S{t + At) - In S(t) = AAt + eaVAl. 



(23) 
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Starting from z , we thus evaluate the expectation value £7[Oi|So] with p = 
2m + l,m 6 N values of z\ centered^ on the mean value E[zi\ = z + AAt 
and which differ from each other of a quantity of the order of ay/At 

z{ = Zq + AAt + jay/At, j = —m, . . . , +m. 

Going on like this, we can evaluate each expectation value E[02\z{] obtained 
from each one of the ziS created above with p values for z<i centered around 
the mean value 

E[z 2 \z{] = z{ + AAt = z + 2 A At + jayfAb. 

Iterating the procedure until the maturity, we create a deterministic grid of 
points such that, at a given time t iy there are (p — l)i + 1 values of Zi, in 
agreement with the request of linear growth. 

This procedure of selection of integration points, together with the calculation 
of the transition probability previously described, is the basis of our path 
integral simulation of the price of a generic option. 



4 Numerical results and discussion 

By applying the results derived in Section 3, we have at disposal an efficient 
path integral algorithm both for the calculation of transition probabilities 
and the evaluation of option prices. In the present section, the application 
of the method to European and American options in the BSM model is illus- 
trated and comparisons with the results obtained with the standard procedures 
known in the literature are shown. 

First, the path integral simulation of the probability distribution of the loga- 
rithm of the stock prices, p(lnS), as a function of the logarithm of the stock 
price, for a BSM-like stochastic model, as given by eq. (2), is shown in Fig. 1. 
The parameters used in the simulation are: So = 100, X = 110, fi = 0.05, 
a — 0.1, t — year and T = 1 year, with 100 time slices. As can be seen, the 
expected lognormal distribution of the stock prices is correctly reproduced by 
the path integral numerical simulation. The plot shows a comparison of the 
calculation of p(lnS) as obtained by means of the two path integral algorithms 
described in Section 3.1. The markers correspond to the naive path integral 
computation of the probability distribution, without variance reduction, for 
10 3 (upper plot) and 10 4 (lower plot) MC iterations. The error bars indicate 

3 Let us recall that between two possible exercise times the probability distribution 
function is gaussian and it is therefore symmetrical with respect to its mean value. 
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Fig. 1. Simulation of the transition probability distribution in the BSM model as a 
function of the logarithm of stock prices via the two path integral methods discussed 
in the text: the naive path integral implementation for 10 3 and 10 4 Monte Carlo 
calls (markers) is compared with the path integral implementation with importance 
sampling (solid line). 

the la - statistical error of the MC calculation. The solid line is the prediction 
for p(lnS) as obtained with the path integral simulation with importance sam- 
pling. In such a case, only two calls for each variables hi are needed to correctly 
fit a gaussian distribution, the numerical error being totally negligible and the 
algorithm very fast, with a typical execution time of a few seconds on a Pentiu- 
mlll 500 Mhz PC. On the contrary, the first path integral implementation is 
much less accurate and CPU time consuming. This is a consequence of the fact 
that, in the path integral simulation with importance sampling, the presence 
of constant integration kernel squeezes to zero the standard estimation error. 
The diagonalization of the tridiagonal matrix M, which is a basic ingredient 
of the efficient path integral algorithm developed, is performed according to 
the standard numerical procedure described in Ref. [16], realized by means of 
the routine F02FAF of the NAG program library [17], while the generation 
of the gaussian variables hi follows from the routine RNORML of the CERN 
program library. 

As previously emphasized, the knowledge of the asset price Si at each interme- 
diate time ti through eq. (22) gives the possibility of applying the procedure 
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Table 1 

Price of an European put option in the BSM model for the parameters t = year, 
T = 0.5 year, r = 0.1, a = 0.4, X = 10, as a function of different stock prices Sq. 
100 time slices are used in the path integral simulation. 



S 


analytical 


binomial tree 


GFDNM 


path integral 


6.0 


3.558 


3.557 


3.557 


3.558 


8.0 


1.918 


1.917 


1.917 


1.918 


10.0 


0.870 


0.866 


0.871 


0.870 


12.0 


0.348 


0.351 


0.349 


0.348 


14.0 


0.128 


0.128 


0.129 


0.128 



of calculation of the transition probability illustrated in Section 3.1 also to 
those types of financial derivatives whose payoff depends on the price path 
of the underlying asset. For such exotic options, the fast computation of the 
transition probability is a basic ingredient of our path integral algorithm as 
an efficient alternative to time-consuming MC simulation. 

Once the transition probability has been computed, the price of an option can 
be computed in a path integral approach as the conditional expectation value 
of a given functional of the stochastic process. For example, the price of an 
European call option will be given by 

+oo 

C = e - r{T - t] J dz f p(z f ,T\zi,t) max[e z f - X,0], (24) 

— oo 

while for an European put it will be 

V = e- r[T - t] J dzfp(z f ,T\zi,t)msx{X - e*',0], (25) 

— oo 

where r is the risk-free interest rate. Therefore just one-dimensional integrals 
need to be evaluated. They can be precisely computed with standard quadra- 
ture rules. In our calculation, the one-dimensional integrals are simply per- 
formed with a standard trapezoidal rule, cross-checked with the routine of 
adaptive integration D01EAF from the NAG library [17]. A sample of the 
results obtained for an European put option in the BSM model is shown 
in Tab. 1. The predictions of our approach, indicated as path integral, are 
compared with results available in the literature, as quoted in Ref. [11]. In 
Tab. 1, the entries correspond to the analytical results, the results by binomial 
trees, and the results of the Green Function Deterministic Numerical Method 
(GFNDM) developed in Ref. [11]. As can be noticed, our results are in per- 
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Solid line — Analytical Markers — Path Integral 




10 20 30 10 20 



Fig. 2. The Greek letters delta (A), gamma (T), vega (V) and theta (i?) for an 
European call option in the BSM model, as functions of the stock price. The lines 
are the analytical predictions, while the markers are the results of our path integral 
algorithm. 

feet agreement with the analytical predictions, while the differences with the 
other numerical procedures are within the 1% level. The errors in our numbers 
as due to numerical integration are not specified being well below the digits 
quoted. 

Further numerical results of our path integral method are shown in Fig. 2, 
where the behaviour of the Greek letters [3] of an European call option in 
the BSM model is plotted as a function of the stock price. The solid lines 
correspond to the analytical predictions, while the markers are the results by 
our path integral approach. The model parameters used in the simulation are: 
X = 10, r = 0.1, a = 0.4, t — year and T = 0.5 year. The four Greeks 
shown are important variables to manage the risk associated to an option and 
are defined as [3] 

a dC t, d2c t r w dC Q 8C 

where C is the price of the European call option. They are computed through 
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Table 2 

Price of an American put option in the BSM model for the parameters t = year, 
T = 0.5 year, r = 0.1, a = 0.4, X = 10, as a function different stock prices So- The 
path integral simulation is performed with 200 time slicings and p = 13 integration 
points. 



S 


finite difference 


binomial tree 


GFDNM 


path integral 


6.0 


4.00 


4.00 


4.00 


4.00 


8.0 


2.095 


2.096 


2.093 


2.095 


10.0 


0.921 


0.920 


0.922 


0.922 


12.0 


0.362 


0.365 


0.364 


0.362 


14.0 


0.132 


0.133 


0.133 


0.132 



standard numerical differentiation of the numerical value C returned by the 
path integral, by means of the NAG routine D04AAF [17]. Also for the Greeks, 
there is perfect agreement between the results obtained via our path integral 
simulation and the analytical predictions. 

To test the reliability of the sampling over the integration points discussed in 
Section 3.2, we present results for the particular case of the price of an Amer- 
ican option in the BSM model in Tab. 2, as obtained with the grid technique 
described in Section 3.2. Comparisons with independent results available in 
the literature are also shown. As can be seen from Tab. 2, there is generally 
a good agreement of our path integral results with those known in the litera- 
ture [11] and obtained by means of the binomial tree, of the finite difference 
method and of the GFDNM method. It is worth noticing that our results 
in the path integral framework require only a few seconds on a Pentiumlll 
500MhZ PC. 



5 The limit of continuum and American options 

In the case of an American option, the possibility of exercise at any time up 
to the expiration date allows to develop, within the path integral formalism, a 
specific algorithm, which, as shown in the following, is precise and very quick. 
Given the time slicing considered in Section 3.2, the case of American options 
requires the limit At — > which, putting a — * 0, leads to a delta-like transition 
probability 

p(z, t + At\z t , t) « 5(z - z t - AM). 
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This means that, apart from volatility effects, the price Zi at time ti will 
have a value remarkably close to the expected value z = + AAt, given 
by the drift growth. Needless to say, if we should substitute the expression 
p(zi,iAt\zi-i, (i — I) At) « S(zi — z) inside the integrals (9), we would neglect 
the role of the volatility and consider only a drift growth of the asset prices. 
In order to take care of the volatility effects, a possible solution is to estimate 
the integral of interest, i.e. 

+00 

EP^] = J dzpizlz^O^), (26) 

—00 

by inserting in eq. (26) the analytical expression for the p(z\zi-i) transition 
probability 

... 1 / (z - Zi-i - AAt) 2 1 

= 7SP exp \ ^At 1 = 

1 f {z-zf \ 

V^At^ P \ 2a 2 At 

together with a Taylor expansion of the kernel function Oi(e z ) = f(z) around 
the expected value z. Hence, up to the second order in z— z, the kernel function 
becomes 

f(z) = f(z) + (z- z)f'(z) + \f'{z\z - zf + 0((z - z) 3 ), (27) 

which, together with eq. (5), yields 

E[O i \S i - 1 ]=m + ^f'(z), + ..., (28) 

since the first derivative does not give contribution to eq. (26), being the 
integral of an odd function over the whole z range. The second derivative can 
be numerically estimated as 

f"(z) = hf(z + 5„) - 2f(z) + f(z - 5 a )], (29) 

with 5 a = O (ay/At), as dictated by the dynamics of the stochastic process. 
It is worth noticing that each expectation value £ , [Oj| J S'j_i] can be now com- 
puted once f(z) = Oi(e Zi - 1+AAt ) and f(z±5 <T ) = Oi(e z *~ 1+AAt±5 °) are known. 
Consequently, if we employ the deterministic grid illustrated in Section 3.2, it 
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Table 3 

Price of an American put option with t = year, T = 0.5 year, r = 0.1, a = 0.4, 
X = 10, as a function different stock prices Sq. The path integral results are obtained 
with S a = 2a y/ At, 300 time slices and p = 3. 



So finite difference binomial tree GFDNM path integral 



6.0 


4.00 


4.00 


4.00 


4.00 


8.0 


2.095 


2.096 


2.093 


2.095 


10.0 


0.921 


0.920 


0.922 


0.922 


12.0 


0.362 


0.365 


0.364 


0.362 


14.0 


0.132 


0.133 


0.133 


0.132 



is enough to put p = 3 to obtain reliable results, provided At is taken suffi- 
ciently small. Actually, the results obtained with this simple semi-analytical 
procedure are shown in Tab. 3, using 5 a = 2a\/At for numerical differentiation 
and 300 time slices. For such path integral evaluation of an American option 
the CPU time is negligible. The results are in nice agreement with those of 
other numerical procedures and in perfect agreement, as we explicitly checked 
for different model parameters, with those quoted in Tab. 2 as obtained with 
the path integral algorithm discussed in Section 3.2. 



6 Concluding remarks and outlook 

In this paper we have shown that the path integral approach to stochastic pro- 
cesses can be successfully applied to the problem of option pricing in financial 
analysis. In particular, an efficient implementation of the path integral method 
has been presented, in order to obtain fast and accurate predictions for a large 
class of financial derivatives, including those with path-dependent and early 
exercise features. The key points of the algorithm are a careful evaluation of 
the transition probability associated to the stochastic model for the time evo- 
lution of the asset prices and a suitable choice of the integration points needed 
to evaluate the quantities of financial interest. Furthermore, a simple and very 
fast procedure to value American options has been derived, by exploiting the 
possibility of continuous exercise up to the expiration date. 

The results of the path integral algorithm have been carefully compared with 
those available in the literature for European and American options in the 
BSM model and found to be in good agreement with the standard numerical 
procedures used in finance. The computational time of the algorithm devel- 
oped is, to the best of our knowledge, competitive with the most efficient 
strategies used in finance. The method is general and it can be quite easily ex- 
tended to cope with other financial derivatives (with path-dependent features) 
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and other models beyond the dynamics of geometric Brownian motion. 

The natural developments of the path integral algorithm here presented con- 
cern the application of the method to value other kinds of quantities of finan- 
cial interest, for which the analytical solution is not available or not accessible, 
and the extension of the method of option pricing to more realistic model of 
the financial dynamics, such as models with stochastic volatility [3,12,18] or 
beyond the BSM Gaussian limit [19-24], in order to search for a better agree- 
ment with the real prices as observed in the real market. A further interesting 
perspective would be using the path integral algorithm as a benchmark to 
train neural networks. 

Such developments are by now under consideration. 
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